namespace il
{

	// Differential Evolution Solver Class
	// Based on algorithms developed by Dr. Rainer Storn & Kenneth Price
	// Written By: Lester E. Godwin
	//             PushCorp, Inc.
	//             Dallas, Texas
	//             972-840-0208 x102
	//             godwin@pushcorp.com
	// Created: 6/8/98
	// Last Modified: 6/8/98
	// Revision: 1.0

	#include <memory.h>
	#if !defined(_DESOLVER_H)
	#define _DESOLVER_H

	#define stBest1Exp			0
	#define stRand1Exp			1
	#define stRandToBest1Exp	2
	#define stBest2Exp			3
	#define stRand2Exp			4
	#define stBest1Bin			5
	#define stRand1Bin			6
	#define stRandToBest1Bin	7
	#define stBest2Bin			8
	#define stRand2Bin			9

	class DESolver;

	typedef void (DESolver::*StrategyFunction)(int);

	class DESolver
	{
	public:
		DESolver(int dim,int popSize);
		~DESolver(void);
	

		// Setup() must be called before solve to set min, max, strategy etc.
		void Setup(double min[],double max[],int deStrategy,
								double diffScale,double crossoverProb);

		// Solve() returns true if EnergyFunction() returns true.
		// Otherwise it runs maxGenerations generations and returns false.
		virtual bool Solve(int maxGenerations);

		// EnergyFunction must be overridden for problem to solve
		// testSolution[] is nDim array for a candidate solution
		// setting bAtSolution = true indicates solution is found
		// and Solve() immediately returns true.
		virtual double EnergyFunction(double testSolution[], bool &bAtSolution) = 0;
	
		int Dimension(void) { return(nDim); }
		int Population(void) { return(nPop); }

		// Call these functions after Solve() to get results.
		double Energy(void) { return(bestEnergy); }
		double *Solution(void) { return(bestSolution); }

		int Generations(void) { return(generations); }

	protected:
		void SelectSamples(int candidate,int *r1,int *r2=0,int *r3=0,
													int *r4=0,int *r5=0);
		double RandomUniform(double min,double max);
		void applyConstraints();

		int nDim;
		int nPop;
		int generations;

		int strategy;
		StrategyFunction calcTrialSolution;
		double scale;
		double probability;

		double trialEnergy;
		double bestEnergy;
	
		double *trialSolution;
		double *lower_bound;
		double *upper_bound;
		double *bestSolution;
		double *popEnergy;
		double *population;

	private:
		void Best1Exp(int candidate);
		void Rand1Exp(int candidate);
		void RandToBest1Exp(int candidate);
		void Best2Exp(int candidate);
		void Rand2Exp(int candidate);
		void Best1Bin(int candidate);
		void Rand1Bin(int candidate);
		void RandToBest1Bin(int candidate);
		void Best2Bin(int candidate);
		void Rand2Bin(int candidate);
	};

	#endif // _DESOLVER_H

	///////////////////////////// Implemetations ///////////////////////////////
	#define Element(a,b,c)  a[b*nDim+c]
	#define RowVector(a,b)  (&a[b*nDim])
	#define CopyVector(a,b) memcpy((a),(b),nDim*sizeof(double))

	DESolver::DESolver(int dim,int popSize) :
						nDim(dim), nPop(popSize),
						generations(0), strategy(stRand1Exp),
						scale(0.7), probability(0.5), bestEnergy(0.0),
						trialSolution(0), bestSolution(0),
						popEnergy(0), population(0)
	{
		trialSolution = new double[nDim];
		bestSolution  = new double[nDim];
		popEnergy	  = new double[nPop];
		population	  = new double[nPop * nDim];

		return;
	}

	DESolver::~DESolver(void)
	{
		if (trialSolution) delete trialSolution;
		if (bestSolution) delete bestSolution;
		if (popEnergy) delete popEnergy;
		if (population) delete population;

		trialSolution = bestSolution = popEnergy = population = 0;
		return;
	}

	void DESolver::Setup(double *min,double *max,
							int deStrategy,double diffScale,double crossoverProb)
	{
		int i;

		strategy	= deStrategy;
		scale		= diffScale;
		probability = crossoverProb;
	
		for (i=0; i < nPop; i++)
		{
			for (int j=0; j < nDim; j++)
				Element(population,i,j) = RandomUniform(min[j],max[j]);

			popEnergy[i] = 1.0E20;
		}

		for (i=0; i < nDim; i++)
			bestSolution[i] = 0.0;

		switch (strategy)
		{
			case stBest1Exp:
				calcTrialSolution = &DESolver::Best1Exp;
				break;

			case stRand1Exp:
				calcTrialSolution = &DESolver::Rand1Exp;
				break;

			case stRandToBest1Exp:
				calcTrialSolution = &DESolver::RandToBest1Exp;
				break;

			case stBest2Exp:
				calcTrialSolution = &DESolver::Best2Exp;
				break;

			case stRand2Exp:
				calcTrialSolution = &DESolver::Rand2Exp;
				break;

			case stBest1Bin:
				calcTrialSolution = &DESolver::Best1Bin;
				break;

			case stRand1Bin:
				calcTrialSolution = &DESolver::Rand1Bin;
				break;

			case stRandToBest1Bin:
				calcTrialSolution = &DESolver::RandToBest1Bin;
				break;

			case stBest2Bin:
				calcTrialSolution = &DESolver::Best2Bin;
				break;

			case stRand2Bin:
				calcTrialSolution = &DESolver::Rand2Bin;
				break;
		}

		return;
	}

	bool DESolver::Solve(int maxGenerations)
	{
		int generation;
		int candidate;
		bool bAtSolution;

		bestEnergy = 1.0E20;
		bAtSolution = false;
	
		for (generation=0;(generation < maxGenerations) && !bAtSolution;generation++)
		{
		
		
			for (candidate=0; candidate < nPop; candidate++)
			{
				(this->*calcTrialSolution)(candidate);
				trialEnergy = EnergyFunction(trialSolution,bAtSolution);

				if (trialEnergy < popEnergy[candidate])
				{
					// New low for this candidate
					popEnergy[candidate] = trialEnergy;
					CopyVector(RowVector(population,candidate),trialSolution);

					// Check if all-time low
					if (trialEnergy < bestEnergy)
					{
						bestEnergy = trialEnergy;
						CopyVector(bestSolution,trialSolution);
					}
				}
				cout << "Generation:"<< generation<<"\n Best Energy: " << bestEnergy << endl;
				cout << "Curr Solution: \n[";
				for(int k = 0 ; k < nDim; ++k)
					cout << bestSolution[k] << ", ";
				cout << "]\n";
			}

		}

		generations = generation;
		return(bAtSolution);
	}

	void DESolver::Best1Exp(int candidate)
	{
		int r1, r2;
		int n;

		SelectSamples(candidate,&r1,&r2);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
		{
			trialSolution[n] = bestSolution[n]
								+ scale * (Element(population,r1,n)
								- Element(population,r2,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::Rand1Exp(int candidate)
	{
		int r1, r2, r3;
		int n;

		SelectSamples(candidate,&r1,&r2,&r3);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
		{
			trialSolution[n] = Element(population,r1,n)
								+ scale * (Element(population,r2,n)
								- Element(population,r3,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::RandToBest1Exp(int candidate)
	{
		int r1, r2;
		int n;

		SelectSamples( candidate, &r1, &r2 );
		n = (int)RandomUniform( 0.0, (double)nDim );

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; ( RandomUniform(0.0,1.0) < probability ) && (i < nDim) ; i++) 
		{
			trialSolution[n] += scale * (bestSolution[n] - trialSolution[n])
								 + scale * (Element(population,r1,n)
								 - Element(population,r2,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::Best2Exp(int candidate)
	{
		int r1, r2, r3, r4;
		int n;

		SelectSamples(candidate,&r1,&r2,&r3,&r4);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
		{
			trialSolution[n] = bestSolution[n] +
								scale * (Element(population,r1,n)
											+ Element(population,r2,n)
											- Element(population,r3,n)
											- Element(population,r4,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::Rand2Exp(int candidate)
	{
		int r1, r2, r3, r4, r5;
		int n;

		SelectSamples(candidate,&r1,&r2,&r3,&r4,&r5);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; (RandomUniform(0.0,1.0) < probability) && (i < nDim); i++) 
		{
			trialSolution[n] = Element(population,r1,n)
								+ scale * (Element(population,r2,n)
											+ Element(population,r3,n)
											- Element(population,r4,n)
											- Element(population,r5,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::Best1Bin(int candidate)
	{
		int r1, r2;
		int n;

		SelectSamples(candidate,&r1,&r2);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; i < nDim; i++) 
		{
			if ((RandomUniform(0.0,1.0) < probability) || (i == (nDim - 1)))
				trialSolution[n] = bestSolution[n]
									+ scale * (Element(population,r1,n)
												- Element(population,r2,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::Rand1Bin(int candidate)
	{
		int r1, r2, r3;
		int n;

		SelectSamples(candidate,&r1,&r2,&r3);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; i < nDim; i++) 
		{
			if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
				trialSolution[n] = Element(population,r1,n)
									+ scale * (Element(population,r2,n)
													- Element(population,r3,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::RandToBest1Bin(int candidate)
	{
		int r1, r2;
		int n;

		SelectSamples(candidate,&r1,&r2);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; i < nDim; i++) 
		{
			if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
				trialSolution[n] += scale * (bestSolution[n] - trialSolution[n])
										+ scale * (Element(population,r1,n)
													- Element(population,r2,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::Best2Bin(int candidate)
	{
		int r1, r2, r3, r4;
		int n;

		SelectSamples(candidate,&r1,&r2,&r3,&r4);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; i < nDim; i++) 
		{
			if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
				trialSolution[n] = bestSolution[n]
									+ scale * (Element(population,r1,n)
												+ Element(population,r2,n)
												- Element(population,r3,n)
												- Element(population,r4,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::Rand2Bin(int candidate)
	{
		int r1, r2, r3, r4, r5;
		int n;

		SelectSamples(candidate,&r1,&r2,&r3,&r4,&r5);
		n = (int)RandomUniform(0.0,(double)nDim);

		CopyVector(trialSolution,RowVector(population,candidate));
		for (int i=0; i < nDim; i++) 
		{
			if ((RandomUniform(0.0,1.0) < probability) || (i  == (nDim - 1)))
				trialSolution[n] = Element(population,r1,n)
									+ scale * (Element(population,r2,n)
												+ Element(population,r3,n)
												- Element(population,r4,n)
												- Element(population,r5,n));
			n = (n + 1) % nDim;
		}

		return;
	}

	void DESolver::SelectSamples(int candidate,int *r1,int *r2,
											int *r3,int *r4,int *r5)
	{
		if (r1)
		{
			do
			{
				*r1 = (int)RandomUniform(0.0,(double)nPop);
			}
			while (*r1 == candidate);
		}

		if (r2)
		{
			do
			{
				*r2 = (int)RandomUniform(0.0,(double)nPop);
			}
			while ((*r2 == candidate) || (*r2 == *r1));
		}

		if (r3)
		{
			do
			{
				*r3 = (int)RandomUniform(0.0,(double)nPop);
			}
			while ((*r3 == candidate) || (*r3 == *r2) || (*r3 == *r1));
		}

		if (r4)
		{
			do
			{
				*r4 = (int)RandomUniform(0.0,(double)nPop);
			}
			while ((*r4 == candidate) || (*r4 == *r3) || (*r4 == *r2) || (*r4 == *r1));
		}

		if (r5)
		{
			do
			{
				*r5 = (int)RandomUniform(0.0,(double)nPop);
			}
			while ((*r5 == candidate) || (*r5 == *r4) || (*r5 == *r3)
														|| (*r5 == *r2) || (*r5 == *r1));
		}

		return;
	}

	/*------Constants for RandomUniform()---------------------------------------*/
	#define SEED 3
	#define IM1 2147483563
	#define IM2 2147483399
	#define AM (1.0/IM1)
	#define IMM1 (IM1-1)
	#define IA1 40014
	#define IA2 40692
	#define IQ1 53668
	#define IQ2 52774
	#define IR1 12211
	#define IR2 3791
	#define NTAB 32
	#define NDIV (1+IMM1/NTAB)
	#define EPS 1.2e-7
	#define RNMX (1.0-EPS)

	double DESolver::RandomUniform(double minValue,double maxValue)
	{
		long j;
		long k;
		static long idum;
		static long idum2=123456789;
		static long iy=0;
		static long iv[NTAB];
		double result;

		if (iy == 0)
			idum = SEED;

		if (idum <= 0)
		{
			if (-idum < 1)
				idum = 1;
			else
				idum = -idum;

			idum2 = idum;

			for (j=NTAB+7; j>=0; j--)
			{
				k = idum / IQ1;
				idum = IA1 * (idum - k*IQ1) - k*IR1;
				if (idum < 0) idum += IM1;
				if (j < NTAB) iv[j] = idum;
			}

			iy = iv[0];
		}

		k = idum / IQ1;
		idum = IA1 * (idum - k*IQ1) - k*IR1;

		if (idum < 0)
			idum += IM1;

		k = idum2 / IQ2;
		idum2 = IA2 * (idum2 - k*IQ2) - k*IR2;

		if (idum2 < 0)
			idum2 += IM2;

		j = iy / NDIV;
		iy = iv[j] - idum2;
		iv[j] = idum;

		if (iy < 1)
			iy += IMM1;

		result = AM * iy;

		if (result > RNMX)
			result = RNMX;

		result = minValue + result * (maxValue - minValue);
		return(result);
	}

	namespace il{
		class DETest1 : public DESolver
		{
		public:
			DETest1(int dim, int pop) : DESolver(dim,_POPULATION_SIZE){}
			double EnergyFunction(double x[], bool &AtSolution)
			{
				//Rastrigin
				double res = 0;
				for(int i = 0 ; i < nDim ; ++i)
					res += x[i] * x[i] - 10 * cos( 2 * CV_PI * x[i] );
				res += 10  *nDim;
				return res;
			}
		};

		class DETest2 : public DESolver
		{
		public:
			DETest2(int dim, int pop) : DESolver(dim,_POPULATION_SIZE){}
			double EnergyFunction(double x[], bool &AtSolution)
			{
				//Rosenbrock
				double res = 0;
				for(int i = 0 ; i < nDim-1 ; ++i)
					res += pow(1-x[i],2) + 100.0 * pow( x[i+1]-pow(x[i],2), 2);
				return res;
			}
		};
	}
}